The Gaussian probability distribution, also called the normal distribution or the bell curve, is an important probability distribution to science, engineering, and the natural world. Its ubiquity is due largely to the Central Limit Theorem, which proves the tendency of many random events such as thermal noise, repeated flips of a coin, or student examination scores (to name a few) to be well described by a Gaussian distribution (see, e.g., A. Papoulis and S. U. Pillai, Probability, Random Variables, and Stochastic Processes, McGraw Hill, Boston, 2002). Gaussian distributions (also referred to as “Gaussians”) can be univariate (one random variable) or multivariate (a vector of jointly Gaussian random variables), and the probability density of a Gaussian at any point can be readily calculated as
                                          p            ⁡                          (              x              )                                =                                    1                                                                    (                                          2                      ⁢                      π                                        )                                                        n                    2                                                  ⁢                                                                          Σ                                                                            1                    2                                                                        ⁢            exp            ⁢                          {                                                -                                      1                    2                                                  ⁢                                                      (                                          x                      -                      μ                                        )                                    T                                ⁢                                                      Σ                                          -                      1                                                        ⁡                                      (                                          x                      -                      μ                                        )                                                              }                                      ,                            (        1        )            where xεRn is a vector with n real valued elements, μεRn is the mean vector defining the center of the Gaussian, and ΣεRn×n is the symmetric, positive semidefinite covariance matrix describing the spread of the Gaussian (in the univariate case, Σ=σ2, the familiar variance term).
In addition to the Central Limit Theorem, the Gaussian receives much attention because of its attractive mathematical properties—for example, unlike many probability distributions, adding and scaling Gaussians yields again a Gaussian. The moment generating function (and characteristic function) of a Gaussian is also a convenient closed-form expression. However, despite its beneficial aspects, the Gaussian presents difficulty in that its cumulative distribution function (cdf) has no closed-form expression and is difficult to calculate. The cdf calculates the probability that a random draw from a Gaussian falls below a certain value (e.g., the probability that a student scores less than a particular value on a test).
In connection with these, cumulative densities represent the probability that a random draw falls in a particular region AεRn. Generally, for a region of interest A in a distribution, where the height of the curve corresponds to the probability density p(x), the cumulative density is the total mass captured above that region A. Mathematically, this is represented as follows
                                          F            ⁡                          (              A              )                                =                                    Prob              ⁢                              {                                  x                  ∈                  A                                }                                      =                                                            ∫                  A                                ⁢                                                      p                    ⁡                                          (                      x                      )                                                        ⁢                                      ⅆ                    x                                                              =                                                ∫                                      l                    1                                                        u                    1                                                  ⁢                                                                  ⁢                                  …                  ⁢                                                                          ⁢                                                            ∫                                              l                        n                                                                    u                        n                                                              ⁢                                                                  p                        ⁡                                                  (                          x                          )                                                                    ⁢                                              ⅆ                                                  x                          n                                                                    ⁢                                                                                          ⁢                      …                      ⁢                                                                                          ⁢                                              ⅆ                                                  x                          1                                                                                                                                                        ,                            (        2        )            where l1, . . . , ln and u1, . . . , un denote the upper and lower bounds of the region A. This cumulative density F(A) generalizes the cdf, as the cdf can be recovered by setting l1= . . . =ln=−∞ (the region A is unbounded to the left). Cumulative densities are important, as is the distribution itself, and the applications are equally widespread, such as those relating to statistics, economics, mathematics, biostatistics, environmental science, computer science, neuroscience and machine learning.
A well-labored body of work has involved investigating methods for calculating these quantities relative to Gaussian cumulative densities. The univariate cdf can be very quickly and accurately calculated using a number of techniques (see, e.g., W. J. Cody, Rational chebyshev approximations for the error function, Math. Comp., pp. 631.638, 1969). These methods are so fast and accurate that the univariate cdf, often denoted ø(•) or a scaled version of the complementary error function erfc(•), is available with machine-level precision (as precise as a digital computer can represent any number) in many statistical computing packages (e.g., normcdf in MATLAB, CDFNORM in SPSS, pnorm in R, to name a few).
However, calculating multivariate Gaussians continues to be challenging. Previous approaches for calculating multivariate cumulative densities have generally required sophisticated numerical integration techniques, which can be undesirable as relative to computational time and power. Importantly, these methods do not produce an analytical form to their answer, so derivatives and other optimizations with respect to cumulative densities cannot be performed. These and other issues remain as a challenge to a variety of methods, devices and systems that use or benefit from Gaussian-based processing.